Under consideration for publication in J. Fluid Mech. 



1 



Universality of shear-banding instability and 
crystallization in sheared granular fluid 

MEHEBOOB ALAM^f, PRIYANKA SHUKLA^ AND STEFAN LUDING^ 

^ Engineering Mechanics Unit, Jawaharlal Nehru Center for Advanced Scientific Research 
Jakkur P.O., Bangalore 560064, India. 
^MultiScale Mechanics, Facuhy of Science and Technology, University of Twente, P.O. Box 217, 7500 AE 

Enschede, The Netherlands 

(Received 15 November 2008) 

The linear stability analysis of an uniform shear flow of granular materials is revisited using 
several cases of a Navier-Stokes' -level constitutive model in which we incorporate the global 
equation of states for pressure and thermal conductivity (which are accurate up-to the maximum 
packing density v„i) and the shear viscosity is allowed to diverge at a density (< v„i), with 
all other transport coefficients diverging at Vm- It is shown that the emergence of shear-banding 
instabilities (for perturbations having no variation along the streamwise direction), that lead to 
shear-band formation along the gradient direction, depends crucially on the choice of the consti- 
tutive model. In the framework of a dense constitutive model that incorporates only collisional 
transport mechanism, it is shown that an accurate global equation of state for pressure or a vis- 
cosity divergence at a lower density or a stronger viscosity divergence (with other transport co- 
efficients being given by respective Enskog values that diverge at Vm) can induce shear-banding 
instabiUties, even though the original dense Enskog model is stable to such shear-banding in- 
stabilities. For any constitutive model, the onset of this shear-banding instability is tied to a 

universal criterion in terms of constitutive relations for viscosity and pressure, and the sheared 
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granular flow evolves toward a state of lower "dynamic" friction, leading to the shear-induced 
band formation, as it cannot sustain increasing dynamic friction with increasing density to stay 
in the homogeneous state. A similar criterion of a lower viscosity or a lower viscous-dissipation 
is responsible for the shear-banding state in many complex fluids. 



1. Introduction 

One challenge in granular flow research is to devise appropriate hydrodynamic/continuum 
models to describe its macroscopic behavior. Rapid granular flows (CampbeU 1990; Goldhirsch 
2003) can be well modeled by an idealized system of smooth hard spheres with inelastic col- 
hsions. The kinetic theory of dense gases has been modified to obtain Navier-Stokes-like (NS) 
hydrodynamic equations, with an additional equation for the fluctuation kinetic energy of par- 
ticles (i.e., granular energy) that incorporates the dissipative-natire of particle collisions. Such 
NS-order hydrodynamic models have been widely used as prototype models to gain insight into 
the "microscopic" understanding of various physical phenomena involved in granular flows. 

The plane Couette flow has served as a prototype model problem to study the rheology (Lun 
et al. 1984; Jenkins & Richman 1985; CampbeU 1990; Sela & Goldhirsch 1998; Alam & Lud- 
ing 2003, 2003a; Tsai, Voth & GoUub 2003; Alam & Luding 2005; Gayen & Alam 2008) and 
dynamics (Hopkins & Louge 1991; McNamara 1993; Tan & Goldhirsch 1997; Alam & Nott 
1997, 1998; Conway & Glasser 2004; Alam et al. 2005; Alam 2005, 2006; Gayen & Alam 
2006; Saitoh & Hayakawa 2007) of granular materials. In the rapid shear flow, the Unear sta- 
bility analyses of plane Couette flow (Alam & Nott 1998; Alam 2006) showed that this flow 
admits different types of stationary and travelling-wave instabilities, leading to pattern forma- 
tion. One such instability is the "shear-banding" instabiUty in which the homogeneous shear flow 
breaks into alternating dense and dilute regions of particles along the gradient direction. This is 
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dubbed "shear-banding" instability since the "nonlinear" saturation of this instability (Alam & 

Nott 1998; Alam et at. 2005; Shukla & Alam 2008) leads to alternate layers of dense and dilute 
particle-bands in which the shear-rate is high/low in dilute/dense regions, respectively, leading to 
"shear-locaUzation" (Varnik etal. 2003). This is reminiscent of shear-band formation in shear- 
cell experiments (Savage & Sayed 1984; Losert et al. 2000; Mueth et al. 2000; Alam & Luding 
2003; Tsai, Voth & GoUub 2003): when a dense granular material is sheared the shearing is 
confined within a few particle-layers (i.e., a shear-band) and the rest of the material remains 
unsheared, leading to the two-phase flows of dense and dilute regions. Previous works (Alam 
& Nott 1998; Alam et al. 2005) showed that the kinetic-theory-based hydrodynamic models are 
able to predict the co-existence of dilute and dense regimes of such shear-banding patterns. 

The above problem has recently been reanalyzed (Khain & Meerson 2006), with reference to 
shear-band formation, with a constitutive model which is Ukely to be valid in the dense limit. 
These authors showed that their 'dense' constitutive model does not admit shear-banding insta- 
bihties of Alam & Nott (1998); however, a single modification that the shear viscosity diverges 
at a density lower than other transport coefficients resulted in the appearance of two-phase-type 
solutions of dilute and dense flows that are reminiscent of shear-banding instabilities. That the 
viscosity diverges stronger/faster than other transport coefficient has also been incorporated pre- 
viously in a constitutive model by Losert et al. (2000) that yields a satisfactory prediction for 
shear-bands in an experimental Couette flow in three-dimensions. 

We revisit this problem to understand the influence of different Navier-Stokes' -order constitu- 
tive models (as detailed in §2) on the shear-banding instabihties in granular plane Couette flow. 
More specifically, we wiU pinpoint how the effects of various models, some of them involving 
the global equation of state and the viscosity divergence (at a density lower than the maximum 
packing), change the shear-banding instabilities predicted by Alam & Nott (1998). For the sake 
of a systematic overview, we will also discuss about the instability results based on special lim- 
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iting cases of these models. One important finding is that a global equation of state {without 
viscosity divergence at a lower density) leads to a shear-banding instability in the framework 
of a "dense" model that incorporates only coUisional transport mechanism. Even with the local 
equation of state, if we use a constitutive relation for viscosity that has a stronger divergence (at 
the same maximum packing density) than other transport coefficients, we recover shear-banding 
instabilities in the framework of the dense-model of Haff (1983). This brings us to a crossroad: 
is there any connection among the results of Alam & Nott, Khain & Meerson, and the present 
work? Is there any universal criterion for the onset of the shear-banding instability in granular 
shear flow? Such an universality for the shear-banding instability indeed exists, solely in terms of 
the constitutive relations, as we show in this paper. Possible connections of the present criterion 
of a lower dynamic friction for the shear-banding state to explain the onset of shear-banding in 
many complex fluids as well as in an elastic hard-sphere fluid are discussed. 

2. Balance equations and constitutive model 

We use a Navier-Stokes-level hydrodynamic model for which we need balance equations for 
mass, momentum and granular temperature: 



Here q = mn = ppU is the mass-density, m the particle mass, n the number density, pp the 
material density and the area/volume fraction of particles; u is the coarse-grained velocity- 
field and T is the granular temperature of the fluid. Note that the granular temperature, T = 
(C^/dim), is defined as the mean-square fluctuation velocity, with C = (c — u) being the 
peculiar velocity of particles and c the instantaneous particle velocity; dim is the dimensionality 




(2.2) 



(2.1) 



(2.3) 
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of the system and here onward we focus on the two-dimensional (dim = 2) system of an inelastic 

"hard-disk" fluid. The flux terms are the stress tensor, P, and the granular heat flux, q; V is the 

rate of dissipation of granular energy per unit volume- for these three terms we need appropriate 

constitutive relations which are detailed below. 

2. 1 . General form of Newtonian constitutive model: Model- A 
The standard Newtonian form of the stress tensor and the Fourier law of heat flux are: 

P = (p-CV-u)I-2mS, (2.4) 
q = -kVT, (2.5) 

where I is the identity tensor and S the deviator of the deformation rate tensor. Here p, /x, C 
and K are pressure, shear viscosity, bulk viscosity and thermal conductivity of the granular fluid, 
respectively. 

Focussing on the nearly elastic limit (e — > 1) of an inelastic hard-disk (of diameter d) fluid, 
the constitutive expressions for p, ji, Q, k and T) are given by 

p{u,T) = Ppfi{u)T, ii{u,T) = ppdf2{u)VT, 
C{i^,T) = ppdfs{u)VT, k{v,T) = ppdf4{u)VT, (2.6) 
D(i.,T) = e^f,{^,e)T^/^, 

where /1-/5 are non-dimensional functions of the particle area fraction i' (Gass 197 1 ; Jenkins & 

Richman 1985): 







(2.7) 
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These constitutive expressions give good predictions for transport coefficients of nearly elastic 
granular fluid up-to a density of w 0.55 (see, figure 2 of Alam & Luding 2003a). In the above 
expressions, xi^) is the contact radial distribution function which is taken to be of the foUowing 
form (Henderson 1975) 

(1 - v/VmY 

that diverges at some finite density v = v^- For the ideal case of point particles (i.e. in one 
dimension, Torquato 1995), we have 1/^ = 1 which is unreahstic for macroscopic grains at very 
high densities; there are two other choices for this diverging density: the random close packing 
density z^r = z^m = 0.82 or the maximum packing density = 7r/2v^ « 0.906 in two- 
dimensions (Torquato 1995); up-to some moderate density {v ^ 0.5), there is no difference in 
the value of x{v) for any choice of = 1 or 0.906 or 0.82. The range of vaUdity of different 
variants of model radial distribution functions is discussed in an upcoming paper (Luding 2008). 

We shall denote the above constitutive model (2.6-2.11), with the contact radial distribution 
function being given by (2.12), as "model- A". Since the stability results do not differ qualitatively 
with either choice of the numerical value for (= 0.82 or 0.906 or 1), we will present all results 
with Vm = 7r/2\/3 in (2.12) (except in figure I2d, see §5.2). 

Now we consider the dilute and dense limits of model-^. It should be noted that each transport 
coefficient has contributions from the 'kinetic' and 'colhsionaF modes of transport: while the 
former is dominant in the Boltzmann Umit (f — > 0), the latter is dominant in the dense hmit {v — > 
Vrn)- For example, the pressure can be decomposed into its kinetic and coUisional contributions: 

p=/+p- = Pp(/fe + /^=)T. (2.13) 

To obtain the constitutive expressions for the dilute and dense regimes, one has to take the ap- 
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propriate limit of all functions /1-/4: 

/l = + /f , /2 = /I + /I, 

(2.14) 

/3 = /I + /I, /4 = /I + /|. 

Based on this decomposition, we have the following two limiting cases of model-^. 

2.1.1. Dilute limit: Model Aq 

For this dilute {u — > 0) model, xi'^) ~^ 1 the constitutive expressions are assumed to 
contain contributions only from the kinetic mode of transport: 

/i = f^{u) = V, h = f^{v) = ^ + ^U, 

/3 = fH'^) = 0, /4 = /Kz.) = ^ + ^>y, (2-15) 

We shall call this "model ^0 "• Note that /s has only the leading-order colUsional contribution, 
since no energy is dissipated in kinetic (free-flow) motion. 

2.1.2. Dense limit: Model Ad 

For this dense model, the constitutive expressions are assumed to contain contributions only 
from coUisional mode of transport: 

/i = /f (i.) = 2u\, h = f^{^) = ^ (1 + f ) Z^'X, 

We shaU call this "model Ad" which is nothing but Haff 's model (1983). 

2.2. Model B: global equation of state (EOS) 

In two-dimensions, a global equation of state for pressure has recently been proposed by Luding 
(2001): 

p/ppT = fi{,y)=iy + iy{P^ + m{iy)[Pdense-P4]), V ^ ly ^ Um = (2.17) 
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where 

Xi\^) - (i-^yi - 128(1-1^)1' 

p _ 2umh:,{y,r.~<^) _ -, (2.18) 

dense — (2y^_j^) ^ ' 

h^{ym-v) = l-0.04(z^„-i^)+3.25(i/™-i/)3, 

m{v) = [1 + exp(-(i^ - t//)/mo)]~^ . 
Here Vf\s the freezing point density, and m{v) is a merging function that selects P4 for v « vf 

and Pdense for V » Vf\ the value of mo is taken to be 0.012 along with a freezing density of 

Vf ~ 0.7. It should be noted that the above functional form of fi{v), (2.17), is a monotonically 

increasing function of v, and it has been verified (Luding 2001, 2002; Garcia-Rojo, Luding & 

Brey 2006) from molecular dynamics simulations of elastic hard-disk systems that the numerical 

values for different constants in ( 12.181 1 are accurate (within much less than 1% except around Vf) 

up-to the maximum packing density. 

Rewriting equation (12.17b as 

h{y)^y + 2v\P{v), (2.19) 
we can define an "effective" contact radial distribution function for pressure: 

X^V) - ^ (^4 + [Pdense - Pa]) , 



Xi + m{v) 



h3.{vm-v) 1 

7, X4 



TT 

with Vm ~ — p. (2.20) 
2v3 



It should be noted here that this functional form of x^(j^) has been verified by molecular dynam- 
ics simulations of plane shear flow of frictional inelastic disks (Volfson, Tsimring & Aranson 
2003). In particular, they found that the simulation data on G{v) — v~)<f{v) agree with the pre- 
dictions of (2.20), however, its divergence appears to occur at the random close packing hmit in 
two-dimensions (;/,„ = j/^ = 0.82). 

Similar to pressure, a global equation for thermal conductivity has been suggested by Garcia- 
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Rojo et al. (2006) which is also accurate up-to the maximum packing density. More specifically, 

the non-dimensional function of density for thermal conductivity (equation (2.10)) is replaced by 

where x'^ {'^) is the "effective" contact radial distribution function for thermal conductivity: 



(2.22) 



_u{l-l//l/m) 2i/ 

with x(i', t'm = 1) being obtained from (2.12) by putting i/^ = 1. 

The model-^ with the above global equation of state for pressure and the global equation 
for thermal conductivity is called "model-i?". Note that the dilute limit of model- i? is the same 
as that of model-j4, however, the dense limits of both models are different in the choice of the 
equation of state and thermal conductivity. 

Now we identify two subsets of model-_B: "model-S^" and "model-B'^" where the former is 
model-i? with a global equation of state for pressure (2. 17) only and the latter is model- -B with a 
global equation for thermal conductivity (2.21) only. As clarified in the previous paragraph, the 
remaining transport coefficients of model- -B are same as in model- A. 

2.3. Model-C: viscosity divergence 

Here we consider a variant of model-A which incorporates another ingredient in the constitutive 
model: the shear viscosity diverges at a density lower than other transport coefficients (Garcia- 
Rojo et al. 2006). This can be incorporated in the corresponding dimensionless function for the 
shear viscosity (2.8): 

/l'M=/2Mfl + :^V (2.23) 



which diverges at a density, v = < i/^, that is lower than the close packing density. The first 
term on the right-hand- side is the standard Enskog term (2.8) and the second term incorporates 
a correction due to the viscosity divergence. For all results shown here, we use =0.71 which 
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was observed for the unsheared case (Garcia-Rojo et al. 2006); note, however, that the precise 
density at which viscosity diverges in a shear flow can be larger than 0.71 (see figure 6 of Alam 
& Luding 2003). 

The model with all transport coefficients as in model-^ but with its viscosity divergence being 
at a lower density is termed as "model-C". Similar to model- A, we can recover the dilute and 
dense Umits of model-C, by separating the kinetic and coUisional contributions to each transport 
coefficient. The dense limit of model-C is, however, reached &\. v = due to the viscosity 
divergence. 

2.4. Model-D: Global equation of state and viscosity divergence 

This is the most general model in which we incorporate: (1) the global equation of state for 
pressure (2.17), (2) the global equation for thermal conductivity (2.21), and (3) the viscosity 
divergence (2.23). The other transport coefficients of model-D are same as in model-^. 

3. Plane shear and linear stability 

Let us consider the plane shear flow of granular materials between two waUs that are separated 
by a distance H: the top wall is moving to the right with a velocity 11^^/2 and the bottom wall 
is moving to the left with the same velocity. We impose no-shp and zero heat-flux boundary 
conditions at both waUs: 



The equations of motion for the steady and fuUy developed shear flow admit the foUowing 
solution: 



n=[u{y),v{y)] = [{Uy,/H)y,% i^{y) = const. =V, T = dP ( ^] ^fp-. (3.2) 



The shear rate, du/dy = Uw/H, is uniform (constant) across the Couette gap, and this solution 
wiU henceforth be called uniform shear solution. Note that if viscosity diverges at = z/^, there 



u = ±Uw/2 and q = at y = ±H/2. 



(3.1) 
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is no uniform shear solution for u > v^, i.e. the uniform shear flow is possible only for densities 

We have non-dimensionalized all quantities by using H, H /Uia and Uw as the reference 
length, time and velocity scales, respectively. The explicit forms of dimensionless balance equa- 
tions as well as the dimensionless transport coefficients are written down in Appendix A. There 
are three dimensionless control parameters to characterize our problem: the scaled Couette gap 
H ^ H/d, the mean density (area fraction) 1^=17 and the restitution coefficient e. Here onwards, 
aU quantities are expressed in dimensionless form. 

3.1. Linear stability 

The stability analysis of the plane shear flow has been thoroughly investigated (Alam & Nott 
1998; Alam 2005, 2006; Alam et al. 2005) using the constitutive models of class A for which 
all transport coefficients diverge at the maximum packing fraction i/ ^ as outlined in §2.1. 
The same analysis is carried out here for a specific type of perturbations that are invariant along 
the streamwise/flow (x) direction, having variations along the gradient (y) direction only. This 
implies that the ^-derivatives of all quantities are set to zero (d/dx{-) = 0) in the governing 
equations. The analysis being identical with that of Alam and Nott (1998), we refer the readers 
to that article for mathematical details. 

Consider the stabiUty of the uniform shear solution ( 13 .21) against perturbations that have spatial 
variations along the y-direction only, e.g. the density field can be written as i^{y, t) — v+v'{y, t), 
with the assumption of small-amplitude perturbations, \iy'{y, t) /v\ << 1, for the linear analysis. 
Linearizing around the uniform shear solution, we obtain a set of partial differential equations: 

with BX=(^,l,^ ■{u',v',T')^{), (3.3) 

where X = {v' , u', v' , T') is the vector of perturbation fields, C is the linear stability operator 
and B denotes the boundary operator (i.e., zero slip, zero penetration and zero heat flux). Assum- 
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ing exponential solutions in time X{y, t) = X{y) exp{u!t), we obtain a differential-eigenvalue 
problem: 

^ A 

^ dy^ ' dy 



ujX = C[ —,—,...] X, with BX^O, (3.4) 



where X{y) — {0,u,v,T){y) is the unknown disturbance vector that depends on y. Here, 
u! = ujr + is the complex frequency whose real part uir denotes the growth/decay rate of 
perturbations and the imaginary part cui is its frequency which characterizes the propagating 
(uii 7^ 0) or stationary (uji ~ 0) nature of the disturbance. The flow is stable or unstable if cj^ < 
or ojr > 0, respectively. 

3.2. Analytical solution: dispersion relation and its roots 

Before presenting numerical stability results (§4), we recall that the above set of ordinary differ- 
ential equations i3.4i admits an analytical solution (Alam & Nott 1998): 

{Hy),T{y)) - (i^i, Ti) cos kniy ± 1/2) (3.5) 
{u{y),v{y)) = (ui,ui)sinfc„(y ± 1/2), (3.6) 

where fc„ = nvr is the 'discrete' wavenumber along y, with n = 1,2,... being the mode number 
that tells us the number of zero-crossing of the density/temperature eigenfunctions along y G 
(—1/2, 1/2). With this, equation (I3.4l i boils down to an algebraic eigenvalue problem: 

AXi^LoXi, (3.7) 

where Xi = {vi, ui,vi, Ti) and ^ is a 4 x 4 matrix. This leads to a fourth-order dispersion 
relation in ui: 

uj'^ + a^uj^ + 02^^ + aiu; + oq = 0, (3.8) 



with coefficients 



0,0 — TjTao4 + -fjwOoQ, Oi — TT3-ai2 + 77^014 + T7B-ai6i 

(3.9) 

O2 — -^0.22 + -3^024, 03 — 130 + -fflOj,2. 
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Here 's are real functions of density (v), temperature (T), and the restitution coefficient (e) 
whose exphcit forms are written down in Appendix B. 

Out of four roots of (13.8b . two roots are real and the other two form a complex-conjugate 
pair It is possible to obtain an approximate analytical solution for these four roots using the 
standard asymptotic expansion for large Couette gaps (H = H/d), with the corresponding small 
parameter being H~^. The real roots have the approximations for large H: 



,(1) 



1 


1 (2) 
ai2 + 022^^ 


(2)^" 
+ 032^^0 









where 

The real and imaginary parts of the conjugate pair, 

^(3.4) ^ ^(3,4) ±i ^(3.4)^ 

have the asymptotic approximations for large H: 



0{H- 



(3.10) 
(3.11) 

(3.12) 

(3.13) 



(3,4) _ 



\ ao4 
IP 

(3,4)^2 ^ _}_ ( ^ 



( 012022. ] 
\ O30 ) 



2ai2 
O {H-^) . 



(3.14) 
(3.15) 



For full models (i.e.. A, B, C and D), it can be verified that ojf'^^ is always negative, making the 
first real root given by (13.101 1. the least-stable mode. However, for the dense models (i.e.. 
Ad, Bd, Cd and Dd), cui^'^^ could be positive, making the travelling waves, given by (13.131 1. the 
least-stable mode at low densities. These predictions have been verified against numerical values 
obtained from spectral method as discussed in §4. 
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4. Stability results: comparison among different models 

For results in this section, the differential eigenvalue problem ( 13.41 ) has been discretized us- 
ing the Chebyshev spectral method (Alam & Nott 1998) and the resultant algebraic eigenvalue 
problem has been solved using the QR-algorithm of the Matlab-software. The degree of the 
Chebyshev polynomial was set to 20 which was found to yield accurate eigenvalues. In princi- 
ple, we could solve (3.8) to obtain eigenvalues, but it provides eigenvalues only for a given mode 
number n = 1, 2, ... in one shot; therefore, one has to solve (3.8) for several n to determine the 
growth rate of the most unstable mode. The advantage of the numerical solution of (13.4b is that 
it provides all leading eigenvalues in one shot. 

4. 1 . Results for model-A and its dilute and dense limits 

As mentioned before, the stability analysis of the uniform shear flow with a 3D-variant (i.e. for 
spheres) of model-A has been performed before (Alam & Nott 1998; Alam 2005, 2006; Alam 
et al. 2005). Even though the results for our 2D-model are similar to those for the 3D-model, we 
show a few representative results for this constitutive model for the sake of a complete, systematic 
study and for comparison with other models; note that the results for model- v4o and modt\-Ad 
are new. 

The phase diagram, separating the zones of stability and instability, in the {v, i?) -plane is 
shown in figure 1 for model-A. The flow is unstable inside the neutral contour, denoted by '0', 
and stable outside; a few positive growth-rate contours are also displayed. For the same parameter 
set, from the respective contours of the frequency, w^, in the [v^ i/) -plane it has been verified that 
these instabilities are stationary, i.e., uji = 0. It is seen that there is a minimum value of the 
Couette gap {H = He ) and a minimum density {v — Vcr) below which the flow remains stable. 
With increasing value of e, the neutral contour shifts towards the right, i.e.. Her increases and 
hence the flow becomes more stable with increasing e. We shall discuss the dependence of e on 
the shear-banding instability and the related instability length scale in §5.2. 




Figure 2. (a) Variation of the growth rate of the least stable mode with H foTv = 0.5 and e = 0.9. (6) 
Density eigenfunctions, !>(j/), of the least stable mode at H = 50, 100 and 150. 
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Figure 3. (a) Phase diagram, showing the positive growth-rate contours, for model-^d (i.e., the dense 
limit of model-A): e = 0.9, Um = 7r/2\/3. (b) Contours of frequency, indicating that the instability in panel 
a is due to travelling waves. 

FigurelUa) shows the variation of the growth rate of the least stable mode, wj. = max w^, with 
Couette gap for i' = 0.5, with other parameters as in figure[T] The kinks on the growth-rate curve 
correspond to crossing of modes n — 1, 2, 3, • • •. This can be verified from figure |2f6) which 
displays density eigenfunctions for three values of Couette gaps H ~ 50, 100 and 150. The 
density eigenfunction at H ~ 50 corresponds to the mode n = 1 (i.e. ~ cos(7r(y ± 1/2)) = 
sin(7ry), see equation (3.5)), the other two at H = 100, 150 correspond to modes n = 2, 3 (i.e., 
i) ~ cos(27ry) and sin(37ry), respectively). In fact, there is an infinite hierarchy of such modes 
as oo which has been discussed before (Alam & Nott 1998; Alam et al. 2005). 

For the dilute limit of model-yl (i.e., model- Aq), there are stationary instabilities at finite 
densities (i/ > 0), and the stability diagram (not shown for brevity) in the (H, i/)-plane looks 
similar to that for model-^ (figure 1), but the range of H over which the flow is unstable is much 
larger. As expected, both model-A and model-Ao predict that the flow is stable in the Boltzmann 
limit (i/ 0). 

Figure 3(a) shows the analogue of figure 1 for model-A^ for which the constitutive relations 
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Figure 4. Variations of (a) the growtii rate and (b) the frequency of the least stable mode with H for 
= 0.1 and e = 0.9. Other parameters as in figure[3] 

are expected to be valid in the dense limit (since the constitutive relations contain the collisional 
part only, §2.1.2); figure 3(6) shows the contours of frequency (corresponding to the least-stable 
mode) in the (H, iy)-plane. The flow is unstable to travelling-wave instabilities inside the neutral 
contour. For this dense model, the crossings of different instability modes (n — 1, 2, • • •) with 
increasing Couette gap H and their frequency variations can be ascertained from figured From 
a comparison between figures 1 and 3, the following differences are noted: 

(1) Model- Ad predicts that the flow is stable in the dense limit which is in contrast to the predic- 
tion of the full model (i.e., model- A) for which the dense flow is unstable. This is a surprising 
result since the kinetic contribution to transport coefficients is small in the dense limit, and hence 
both model-A and model- A^ are expected to yield similar results. 

(2) There is a travelling-wave (TW) instability at low densities (v < 0.3) for model-A^ which 
is absent in model-A. Since model- is devoid of the kinetic modes of momentum transfer 
and hence not applicable at low densities, we call the TW-modes in figure 3(a) as "anomalous" 
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Figure 5. Phase diagram, showing the positive growth-rate contours, for model-i?'^ with a global 



modes and they vanish when both kinetic and coUisional effects are incorporated as in the fuU 
model- A. 

One conclusion that can be drawn from the results of three variants of model-^ is that the 
choice of the constitutive model is crucial for the prediction of shear-banding instabiUty. We 
shall come back to discuss this point in § 5. 



Figure 5 shows a phase-diagram in the [H, z/)-plane for model B** ; the flow is unstable inside the 
neutral contour and stable outside. Recall that this model is the same as model-^ for < Vf, with 
the only difference being that we use a global equation for thermal conductivity which is valid 
up-to the maximum packing density (equation 2.21). This instabihty is stationary and the other 
features of stabihty diagrams remain the same (as those in figure 1 for the standard model- A), but 
there is a dip on the neutral contour at the freezing density v = vj. For its dense counterpart, the 
model- does not predict any instabihty at large densities but has travelhng-wave instabihty at 



equation for thermal conductivity: e = 0.9, Vm = 7r/2v'3, Vf = 0.7. 



4.2. Results for model- B: influence of global equation of state 
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H H 

Figure 6. Phase diagram, showing the positive growth-rate contours, for model-S*' with a global 
equation for pressure: e = 0.9, fm = 7r/2v^, vf = 0.7. (a) Full model; (6) dense limit. 

low densities (the corresponding stability diagram looks similar to that in figure 3(a) for model- 
Ad and hence is not shown). 

When a global equation of state for pressure is incorporated (i.e., model- B^", see §2.2), the 
phase-diagram in the (H, z/)-plane looks markedly different, especially at > Pf, as seen in 
figures 6(a) and 6(6). (Recall that the model-S^ is same as model- A, except that we use a global 
equation of state for pressure, equation (2.17).) In figure 6(b), the neutral stabiUty contour con- 
tains a kink at i/ « 0.3 and there are two instabihty-lobes: the lower instability lobe is due to 
travelling-waves and the upper-one is due to stationary-waves. (Similar to model- A^, the low 
density TW instability in figure 6(b) is dubbed "anomalous" since model- is not vaUd at low 
densities.) It is interesting to note in figure 6(6) that for the dense hmit of model- (i.e., model- 
B^) the flow remains unstable to the "stationary" shear-banding instabiUty up-to the maximum 
packing density. This observation is in contrast to the predictions of model- (figure 3a) and 
model- B^. Therefore, we conclude that within the framework of a dense model the global equa- 
tion of state for pressure induces shear-banding instabilities at large densities. 

When both the global equations for pressure and thermal conductivity are incorporated, the 
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Figure 7. Phase diagram for model-C with viscosity divergence: e = 0.9, Vm = 7r/2V^, = 0.71. 

(a) Full model; (6) dense limit. 

phase diagrams in the {H, z/)-plane look qualitatively similar (not shown) to those for model- 
s'" as in figure 6, with the only difference being slightly higher growth rates for the least-stable 
mode. It is noteworthy that with the global EOS, the flow becomes unstable to shear-banding 
instabilities at very small values of H for v ^ Uf. 

From this section, we can conclude that within the framework of a "dense" constitutive model 
(that incorporates only coUisional contributions to transport coefficients, §2.1.2), a simple mod- 
ification with a global equation of state for pressure induces new shear-banding instabihties at 
large densities; however, a similar modification with a global equation of state for thermal con- 
ductivity does not induce any new instability. 

4.3. Results for model-C and model-D: influence of viscosity divergence 

As discussed in §2.3, model-C is the same as model-^, with the viscosity divergence being at 
a lower density u = < Um- On the other hand, model-D is the most general model that 
incorporates the viscosity divergence at a lower density i' = I'fj. < i^m (equation 2.23) as well 
as global equations for pressure (equation 2.17) and thermal conductivity (equation 2.21). The 
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stability results for these two models are found to be similar, and hence we present results only 
for model-C. 

Figures 7(a) and 1(b) show phase-diagrams in the (H, z/)-plane for model-C and its dense vari- 
ant model-Cd, respectively; the results are shown up-to the viscosity divergence since uniform 
shear is not a solution for v > v^. For the full model in figure 7(a), the phase-diagram looks 
similar to those for model-A and model- i?. For its dense counterpart in figure 7(6), the phase- 
diagram is similar to that for model-B^ and model-B^ in the sense that all three models support 
shear-banding instabilities at large densities. Therefore, for model-Cd, the viscosity divergence 
induces shear-banding instabilities at large densities. This prediction is in tune with the results of 
Khain & Meerson (2006) whose model is similar to our model-C^ (see § 5.3). 

Within the framework of a "dense" model (§2.1 .2), therefore, we can conclude about the emer- 
gence of shear-banding instabilities at large densities: (1) a global equation of state for pressure 
alone can induce shear-banding instabilities; (2) a viscosity-divergence at a density, v = v^, 
lower than the maximum packing density alone can induce shear-banding instabilities. 

5. Discussion: Universality of shear-banding instability and crystallization 

5.1. An universal criterion for shear-banding instability 

Since the shear-banding instability is a stationary (uji = 0) mode, this instability is given by one 
of the real roots (equation (3.10)) of the dispersion relation. The condition for neutral stabiUty 
(ujr — 0) can be obtained by setting uj — in the dispersion relation ( I3.8l i: 



ao = kl/H'^^^, 



(5.1) 



where 
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For an instability to occur at a given density, there must be a range of 'positive' discrete wave- 
numbers, i.e., kn/H > 0, which is equivalent to 

*2 > 0, 

since is always positive. The expression for can be rearranged to yield (Alam 2006) 

> 0, (with /i^ > 0) (5.2) 



V fl 

which must be satisfied for the onset of instability. (It should be noted that (15 .21) provides a 
necessary condition for instability, but the sufficient condition is tied to the thermal-diffusive 
mechanism that leads to an instabiUty length scale (15.11 ) which is discussed in §5.2.) The term 
within the bracket in ( 15.2b is the ratio between the shear stress, P^y = fij, and the pressure, p, 
for the plane Couette flow: 

. _ P.y ^ ) ^ f2VT ^ /2 _ /2 _ 

where 7 = du/dy is the local shear rate (which is a constant for uniform shear flow). This 
is nothing but the dynamic friction coefficient of the shear flow, which must increase with in- 
creasing density for the shear-banding instability to occur Note that as per Navier-Stokes-level 
description, the steady fully developed plane Couette flow admits solutions for which the shear 
stress and pressure are constants across the Couette gap. Hence, the dynamic friction coeffi- 
cient, 13d = /i(du/dj/)/p — /17/p, is a /^oi/f/on-mc/e/^enc/enf order-parameter for both "uniform" 
(7 = const) and "non-uniform" (7 = 7(2;)) shear flows. 

For model- A, the variation of (3d with v is non-monotonic as shown by the solid line in figure 8: 
in the dilute limit (3d is large whose value decreases with increasing density till a critical density 
V = Vcr is reached beyond which (3d increases. Recall that the onset of shear-banding instability 
corresponds to > 0, denoted by the square-symbol in the inset of figure 8. It is clear from 
this inset that the full model is unstable to shear-banding instability for v > Vcr, the dilute model 
is unstable at any density (since ^ > for > 0), but the dense model is stable for aU densities 
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Figure 8. Variation of the "dynamic" friction coefficient, Pa = m/p, with density for model-j4: e = 0.9, 
Vm = 7r/2^/3. The inset shows the variation of ^ with v: the onset of shear-banding instability corre- 
sponds to > 0, denoted by the square-symbol. 

(since ^ = 0). It is worth pointing out that, for the full model- A, the critical density for the 
onset of the shear-banding instability (cf. figure 1) is Vcr ~ 0.31733 which is independent of 
the restitution coefficient e as confirmed in figure 9. The lower branch of the neutral contour of 
figure 1 asymptotically approaches this critical density as — > oo. 

For model-B, the variation of /3d with density {v) is shown in figure 10, and the inset shows 
the variation of dPa/dv with i^. (Recall that this model incorporates the global equations of 
states for pressure and thermal conductivity.) In contrast to the 'stable' model- A^, the model- 
ed (i.e. the dense variant of model-B) is unstable to shear-banding instabilities for all densities 
(figure 66) since ^ > 0. For the full model-S, the critical density for the onset of shear-banding 
instability (cf. figure 6a) is fcr ~ 0.2351 which is also independent of the restitution coefficient 
e (as in figure 9 for model- A). For models C and D, this critical density is ~ 0.2849 and 
for ~ 0.2207, respectively. 

It should noted that there is a constraint {fii, > 0) with our instability criterion (5.2), i.e.. 
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Figure 9. The effect of restitution coefficient e on the variation of with v for model A. 




Figure 10. Same as figure[8] but for model B. 
fi should be a monotonically increasing function of density for the instabiHty to occur. This 
constraint on /i, see equations (2.7) and (2.19), is satisfied for all four models. The consequence 
of a possible non-monotonic /i is that the shear flow remains stable over a small density range 
between freezing and melting. However, the increasing dynamic friction with increasing density 
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(5.2) still remains the criterion for the onset of the shear-banding instabiUty. Even though for 

an unsheared elastic system, /i is non-monotonic (Luding 2001) between freezing and melting 

density, we retain equation (2.19-2.20) for /i since this functional form has been shown to hold 

up-to the random close-packing density in simulations of plane shear flow of inelastic hard-disks 

(Volfsonera/. 2003). 

5.1.1. Shear-banding state and crystallization 

For all four models (A, B, C and D), the shear flow is stable and the system remains ho- 
mogeneous at low densities, but becomes unstable to shear-banding instabihty beyond a critical 
density {u > Ucr)- This is due to the fact that for v > u^r the flow cannot sustain the increas- 
ing dynamic friction {(id) and hence breaks into alternating layers of dilute and dense regions 
along the gradient direction. For v > Vcr, the associated "finite-amplitude" bifurcated solution 
corresponds to a lower shear stress, or, equivalently, a lower dynamic friction coefficient. This 
has been verified (Alam 2008) numerically by tracking the bifurcated solutions of the associated 
steady nonUnear equations. 

A representative set of such nonlinear bifurcated solutions for the profiles of density 
granular temperature T{y) and stream-wise velocity u{y) are displayed in figure 11 for three 
values of the Couette gap H = 50, 75 and 125 for model-A; the related numerical procedure is 
the same as described in Alam et al. (2005). For this plot, the mean density is set to 0.5 and the 
restitution coefficient is 0.9; the corresponding growth-rate variation of the least stable mode can 
be ascertained from figure 2(a). These nonlinear solutions bifurcate from the n = 1 mode (see, 
equations 3.5 and 3.6) of the corresponding linear stability equations, and there is a pair of non- 
linear solutions for each H due to the symmetry of the plane Couette flow (Alam & Nott 1998). 
For mode n = 1, the density is maximum at either of the two walls, and this density-maximum 
approaches the maximum packing density (f^ w 0.906) at H = 125 (soUd line in figure 11a) 
for which we have the coexistence of a "crystaUine" zone and a dilute zone, representing a state 
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of phase separation. Within the crystalline zone, the granular temperature approaches zero (fig- 
ure lib) and so does the shear rate (figure 11c). It is noteworthy that the shear-rate is almost 
uniform and localized within the dilute zone. The resulting two-phase solution is called a "shear- 
band" (or, shear-localization) since the shearing is confined within a band of an agitated dilute 
region that coexists with a denser region with negligible shearing (i.e., a crystalhne-region at 
large enough Couette gap). 

At -ff = 50 with parameters as in figure 11, there are two possible solutions: a "uniform- 
shear" state, with the "dynamic" friction coefficient being (3d ~ 0.26965, which is unstable, and 
one of two "shear-band" solutions for which (3^ w 0.2672. The selection of the stable branch is 
determined by the value of the dynamic friction coefficient being the lowest among all possible 
solutions, and therefore the equiUbrium state of the flow (at H = 50) corresponds to the shear- 
banding state of a "lower" dynamic friction (Alam 2008). 

For higher-order modes (n = 2, 3, • • •), the shape of the nonlinear density/temperature/velocity 
profiles can be ascertained from (3.5) and (3.6). For example, the density profiles for modes 
n = 2 and 3 would look like the corresponding density eigenfunctions in figure 2(6). In fact, the 
solution for the first mode (n = 1) serves as a "building-block" of solutions for higher modes 
which has been clarified previously (Alam & Nott 1998; Nott et al. 1999; Alam et al. 2005). 

At this point, we can say that an accurate constitutive model over the whole range of densities 
(that incorporates both kinetic and colhsional modes of transport mechanisms) should be used 
since the dilute and dense regimes can coexist even at a moderately low mean density. Details 
of such shear-banding solutions for other models {B, C and D) and their stabiUty as weU as the 
related results from particle-dynamics simulations will be considered in a separate paper. Such 
an exercise will help to identify and tune the best among the four models. 
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Figure 1 1. Nonlinear shear-banding solutions for model-A for mode n = 1 with v — 0.5 and e — 0.9: 
(a) density, (b) granular temperature and (c) stream-wise velocity. 

5.2. Instability length scale and the ejfect of dissipation 

While our instability criterion, given by equation ( |5.2| l, yields a critical value for the density 
above which the uniform shear flow is unstable to the shear-banding instability, it does not say 
anything about the related instability length scale below which the flow is stable (cf. figures 1-7). 
This issue of a dominant instability length scale is tied to the underlying diffusive mechanisms in 
a granular fluid, offered by the pseudo-thermal conductivity in the energy balance equation (2.3). 
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We have mentioned in §4 that the neutral contour in the (H, i/)-plane (such as in figures 1-7) 
shifts towards right with increasing restitution coefficient (e), and hence the flow becomes more 
stable in the elastic limit. In fact, the dependence of the neutral contour on e can be removed if 
we define a normalized Couette gap as 



which can be thought of as an "instability length scale" (Tan & Goldhirsch 1997; Alam & Nott 
1998). This length scale appears directly from an analysis of the equation for the neutral contour 



where /5(z^, e) = (1 — e^)/5o(!^), and /4(j^) is related to pseudo-thermal conductivity as in (2.6). 
This specific functional dependence of the "instability length scale" on the restitution coefficient 
is also due to the dependence of the thermal conductivity on the granular temperature which 
implicitly depends on e (equation (3.2)): T ~ f2/ ~ (1 — e^)~^. 

In terms of the above instability length scale (|5.4| l. the renormalized stability diagrams in the 
(H*, j/)-plane are displayed in figures 12(a), 12(6), 12(c) for model A, B, C and D, respectively. 
In each panel, we have superimposed three neutral contours for e — 0.9, 0.99 and 0.999 which are 
indistinguishable from each other due to the underlying scaUng ( 15.5b with e. In figure 12(6), we 
compared the neutral contour of the full model-B with those for model-iJ^" (that incorporates the 
EOS for pressure) and model-i?" (that incorporates the global EOS for thermal conductivity). A 
global equation for thermal conductivity has little influence on the stability diagram, but a global 
equation for pressure significantly enlarges the domain of instability in the {H* , i^)-plane. 

Comparing figure 12(a) with figure 12(6), we find that there is a significant difference be- 
tween the predictions of model-yl and model- i?, especially in the dense limit. In particular, with 



H* = H^/l - e2 



(5.4) 



(EB: 



H^^i = = hH^lh = /5o(l - e2)i/V/4 




(5.5) 
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Figure 12. Renormalized stability diagrams in the {H*, z/)-plane, showing the neutral contour that sepa- 
rates stable and unstable regions: (o) model A; (b) model B; (c) model C and D; (d) model A with different 
x(j') as explained at the end of §5.2. The outermost dotted curve in panel (d) is discussed at the end of §5.3. 



In each panel, the abscissa has been renormalized as H* = Hy'l — e^. Note different range of vertical 
axis in each panel. 

accurate equations of state as in model-B, the dense flow is unstable to shear-banding instabiUty 
even for small values of the Couette gap. This is an important issue beyond the square packing 
density v > 7r/4 (in two-dimensions) for which the flow must reorganize internally such that a 
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part of the flow forms a layered crystalline structure or banding (Alam & Luding 2003), thereby 
aUowing the material to shear. This reconciles weU with the predictions of model-S, but not 
with model-^ (which uses the standard Enskog expressions for all transport coefficients, §2.1) 
for H* < 30 (figure 12a); therefore, using an accurate equation of state over the whole range of 
densities should give reasonable predictions for shear-banding solutions. 

In this regard, the predictions of model-C (upper curve in figure 12c) and model-D (lower 
curve in figure 12c) are also consistent for dense flows since the shear-banding instabilities persist 
at smaU values of H for these models too. It may be recalled that these two models (C and D), 
with viscosity divergence at = < I'm, do not admit "uniform" shear solution at large 
densities u > v^. However, the related shear-banding solutions at < i^^ can be continued to 
higher densities v > u^hy embedding the present problem into the uniform shear case such that 
the shear work is vanishingly small (Khain & Meerson 2006). 

Following one referee's suggestion, we briefly discuss possible effects of Torquato's (1995) 
formula for the contact radial distribution function: 

= ^11^. for Q^.^.f, ^^^^ 

= '\T!;;?r>^^\ for u.^.^u,, 
which is known to be vahd for an elastic hard-disk fluid over a range of densities up-to the random 
packing limit Vr = 0.82. When (5.6) is used instead of (2.12) in model- A, the neutral stability 
curve in the {H* , j')-plane follows the thick line in figure 12(d). For the sake of comparison, we 
have also superimposed the neutral stability curve of model-^, denoted by the thin line, with 
x(z^) being given by (2.12) with v^ = v^. Clearly, a larger range in the dense regime is unstable 
with Torquato's formula (5.6), but the overall instability characteristics (the stationary nature of 
instability, the magnitude of growth rate, etc.) remain similar for both (2.12) and (5.6). We have 
also verified that the stability diagram looks similar to that in figure 5 (except for the presence of 
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a discontinuity on the neutral contour at = Vf, resulting in two instability lobes) when (2.20) 

is used as an effective contact radial distribution function for all transport coefficients. 

5.3. Discussion of some "ultra- dense" constitutive models 

Focussing on the "ultra-dense" regime with volume fractions close to the close-packing density 
(i/ — > v,n), the dense limit of model- A (i.e. model in § 2. 1 .2) can be simplified by replacing v 
with Vm and retaining the dependence of the /fs on v via the corresponding dependence of the 
pair correlation function. For this ultra-dense regime the constitutive expressions are: 

This model is devoid of shear-banding instabihties since it can be verified that 

— — =0, V V > Q. 
av 

This is similar to the predictions of model-A^ (see the dot-dash line in figure 8). 

The constitutive model of Khain & Meerson (2006) can be obtained from (15.71 1 by replacing 
the contact radial distribution function for viscosity by 

„ , , / 0.037 \ , , 

that diverges aiv = as in (2.23); the exact density at which viscosity diverges is not important 
here. Interestingly, Khain & Meerson found two-phase-type solutions using their model. Since 
viscosity diverges at a lower density (and hence "faster") than other transport coefficients, it is 
straightforward to verify that our instabihty criterion, 

— — > 0, V V > Vcr, 

av 
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holds for this model. Recently, Khain (2007) also found two-phase-type solutions using a con- 
stitutive model (a modified model of Khain & Meerson 2006) which is similar to our model- _D, 
and the predictions of his model agree well with particle simulation results; for this modified 
model too, our instability criterion i5.2i holds. Therefore, the two-phase-type solutions of Khain 
& Meerson (2006, 2007) are directly tied to the shear-banding instabilities of Alam et al. (1998, 
2005, 2006) via the universal instability criterion (I5.2l l. More specifically, both belong to the 
same class of constitutive instability (Alam 2006). 

The constitutive model of [Losert et al. (200())| can be obtained from ( I5.7l i by using the follow- 
ing functional form for the contact radial distribution function 

X = (1 - I'/l^r)'^ , 

that diverges at the random close packing limit (I'r), for all transport coefficients except the shear 
viscosity, fi, that has a "stronger" rate of increase near i/r'. 

This choice of viscosity satisfies our instability criterion, d/3d/di^ > 0, and, therefore, the model 

of Losert et a/. would yield shear-banding type solutions which is again tied to the increase of 

"dynamic" friction with density for the uniform shear state. 

To illustrate the quantitative effect of a stronger viscosity divergence on the shear-banding 

instability, the neutral stability contour for model-y4 (with a stronger viscosity divergence, see 

below) is shown in figure 12(d), denoted by the outermost dotted curve. For this case, the con- 

stititutive model is the same as the full model-A (i.e. all transport coefficients diverge at the 

same density i^,,,) but its viscosity function f2{i^) in (2.8) is calculated using a radial distribution 

function that has a stronger divergence than all other transport coefficients: 

^, . ^ 1 - 7iy/m 
(1 - iy/iy„i)i 

with its exponent g = 2.25 > 2. (It should be noted that this specific functional form of x'^ may 
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not be correct quantitatively at all densities, but it has simply been chosen to illustrate the possible 

effects of a stronger viscosity divergence on instabilities.) Comparing the dotted contour in figure 

12(d) with the thin-solid contour for model- A, we find that a larger range in the (v, i7*)-plane 

is unstable for the case of a stronger viscosity divergence. With further increase of q, the neutral 

contour shifts towards the left to cover smaller values of H, leading to even larger instability 

region in the (ly, i/*)-plane. In either case, however, the nature of the shear-banding instability 

remains the same and we do not find any new instability as emphasized before. 

5.4. Limit of elastic hard-sphere fluid: dissipation versus effective shear rate 

Naively extrapolating the instability length-scale ( 15. 5t to the elastic limit (e 1) of atomistic 
fluids results into H ^ oo that corresponds to an infinite system for shear-banding to occur 
in an atomistic fluid. This is in contrast to molecular dynamics simulations of sheared "elastic" 
hard-sphere fluid (Erpenbeck 1984) for which a shear-induced ordering phenomenon has been 
observed at moderate densities (much below the freezing density). Similar observations of such 
banding have been made in simulations for continuous potentials too (see, Evans & Morriss 
1986). The key to resolve this apparent anomaly lies with the fact that the elastic limit (e = 1) 
is singular since the collisional dissipation vanishes. To achieve a steady-state in simulations 
of a sheared atomistic fluid, thermostats are used to take away energy from the system that 
compensates the production of energy due to shear-work (P : Vu). Otherwise the system would 
continually heat up, leading to an infinite temperature. Hence, the collisional dissipation in a 
granular fluid can be seen to play the role of a thermostat in an atomistic fluid. 

Equating the dissipation term (either due to a thermostat in an elastic fluid, or due to inelastic 
collisions in a granular fluid) with the shear-work, we obtain a scaUng relation for temperature 
with the shear rate and the restitution coefficient: 
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where T is the dimensional temperature, and 

defined as an "effective" shear rate. For an elastic fluid, this effective shear rate, 7*, is used to 
normalize the temperature, and hence a similar criterion ( 15.21 ) is likely to hold for the onset of 
shear-banding instability in an elastic fluid too. The dependence of the effective shear rate ( 15.91 ) 
with inelasticity suggests that the shear-banding in atomistic fluids is likely to occur at large shear 
rates, a prediction that agrees with Erpenbeck's (1984) simulations. This needs to be checked by 
determining the analytical expression for the thermostat term (which might depend on the choice 
of the thermostat) in the energy equation. 

While the Erpenbeck's ordering transition has been explained (Kirkpatrik and Nieuwoudt 
1986; Lutsko & Dufty 1986) as an instability of the "unbounded" shear flow of an elastic fluid, 
using Navier-Stokes-level equations with wave-vector-dependent transport coefficients (i.e. gen- 
eralized hydrodynamics), the latter work by Lee et al. (1996) has identified a long-wave insta- 
bility (with perturbations along the gradient direction only) in the uniform shear flow for all 
densities. In particular, Lee et al. showed that the Navier-Stokes-level constitutive model is the 
"minimal" model to predict the robustness of this instability. The possible connection of this 
instability with the present work needs to be investigated in the future. 

5.4. 1 . Shearbanding criterion in a molecular fluid: Loose and Hess ( 1989) 

We close our discussion by recalling a similar instability criterion for an "ordering" transi- 
tion in a dense molecular fluid (Loose and Hess 1989) and its connection (Alam 2006) with 
our instability criterion i5.2\ . along with a more general criterion for shear-banding in a shear 
thinning/thickening fluid. 

Using a non-Newtonian constitutive model. Loose and Hess (1989) have derived a criterion 
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for the onset of shear-banding in a dense molecular fluid: 

/ dpyA ( dpyy\ ( dpyA ( Opyy 



\ J \ dv J ^ \ diy ) \ ) ' (5 0) 

where pyy and py^ are the normal and shear stresses, respectively. Assuming the following func- 
tional dependence of pyy and pyx with density {v) and shear rate (7), 

Pyy =Pyyi^)fyvil) and Pyx ^ plA^)fyx{l), (5.11) 
the above instability criterion simplifies to 

^fyx \ ( p '^Pyy \ ^ ( f '^Pyx \ f n '^fyy 
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For a granular fluid, the shear-rate dependence of stresses follows the well-known Bagnold 
scaling: 

fyyil)-Tr^j^ and fyxh) - iVT ^ (5.13) 

where we have used the relation of the granular temperature with the shear rate, T ~ 7^. Substi- 
tuting ( 15.13b into ( 15.121 ). the Loose-Hess instability criterion boils down to (Alam 2006) 



d / Pyx 

dr. Uo 



^ 0. (5.14) 



The term within the bracket is the dynamic friction coefficient of a fluid, and hence the onset of 
instability is again tied to the increasing value of this dynamic friction coefficient with increasing 
density. This is same as our shear-banding instabihty criterion (15. 3t . For a more general case, the 
shear-rate dependence of stresses can be postulated as 

fyyh) = 7'" 

(5.15) 

fyxil) = 7"+\ 

where the index n is a measure of shear-thickening (n > 0) or shear-thinning (n < 0) behaviour 
of the fluid. With this, the shear-banding instability criterion boils down to 



dp^ (5 16) 
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6. Conclusion and Outlook 

To conclude, we showed that by just knowing the constitutive expressions for pressure and 
shear viscosity, one can determine whether any Navier-Stokes' -level constitutive model would 
lead to a shear-banding instability in granular plane Couette flow. The onset of this stationary 
instability is tied to the increasing value of the "dynamic" friction coefficient, /3d — fij/p (where 
/i, p and 7 = du/dy are the shear viscosity, pressure and shear rate, respectively), with increasing 
density for ly > Vcr (equation (15.21) ): the "homogeneous" shear flow breaks into alternating layers 
of dilute and dense regions along the gradient direction since the flow cannot accommodate the 
increasing friction to stay in the homogeneous state. For v > Vcr, the associated "nonlinear" 
shear-band solution corresponds to a lower shear stress, or, equivalently, a lower dynamic friction 
coefficient (Alam 2008). In other words, the sheared granular flow evolves toward a state of 
"lower" dynamic friction, leading to shear-induced band formation along the gradient direction. 
Note that the dynamic friction coefficient, j3d ~ 1^1 /p, is a position-independent order-parameter 
for both "uniform" (7 — const) and "non-uniform" (7 = 7(1/)) shear flows. 

In the framework of a "dense" constitutive model that incorporates only collisional transport 
mechanism (i.e. Haff's model, 1983), we showed that an accurate global equation of state for 
pressure or a viscosity divergence at a lower density (with other transport coefficients being given 
by respective Enskog values) can induce shear-banding instabilities, even though the original 
dense Enskog model is stable to such shear-banding instabilities. Since the prediction of the 
shear-banding instability depends crucially on the form of the constitutive relations, we need to 
use accurate forms of constitutive expressions over the whole range of density that incorporate 
both kinetic and collisional transport mechanisms. The latter statement is important since the 
dilute and dense regimes coexist even at a low mean density when the uniform shear flow is 
unstable to shear-banding instability. The resulting nonlinear shear-banding solutions of all four 
models {A, B, C and D) and their stability as well as the related results from particle-dynamics 
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simulations will be considered in a separate paper This will help to identify the best among these 

four constitutive models, or will point towards new models. In this regard, it is recommended that 

the particle-dynamics simulations be used to find out accurate expressions (vahd over the whole 

range of density) for all transport coefficients. 

We established that the two-phase-type solutions of Khain & Meerson (2006) are directly 
related to the shear-banding instabilities of Alam et al. (1998, 2005, 2006) via the universal 
instabiUty criterion i5.2i . and both belong to the same class of constitutive instability (Alam 
2006). In particular, the instabilities arising out of non-monotonicities of constitutive relations 
with mean-fields (e.g. the coil-stretch transition is tied to non-monotonic stress-strain curve; see, 
de Gennes 1974) are of constitutive origin and hence dubbed constitutive instability. The same 
universal criterion (15.2b also holds for the constitutive model of Losert et al. (2000), thereby 
yielding such two-phase-type solutions in their model of plane shear flow. 

The onset of the ordering transition of Erpenbeck (1984) in a sheared "elastic" hard-sphere 
fluid (which is close to our granular system) is accompanied by a decrease in viscosity and hence 
a lower viscous dissipation. Therefore, similar to the sheared granular fluid, the state of lower 
viscosity/friction is the preferred equilibrium state for a sheared atomistic fluid. Our instability 
criterion ( 15.2b seems to provide a unified description for the shear-banding phenomena for the 
singular case of hard-sphere fluids if we relate the collisional dissipation to a thermostat, leading 
to an "effective" shear rate. This possible connection needs to be investigated further from the 
viewpoint of a constitutive instability of the underlying field equations. 

The shear-banding phenomenon is ubiquitous in a variety of complex fluids under non-equilibrium 
conditions: wormlike micelles (Spenley, Gates & McLeish 1993), colloidal suspensions (Hoff- 
man 1972; Ackerson & Glark 1984), model glassy material (Varnik et al. 2003), suspensions of 
rod-like viruses (Lettinga & Dhont 2004), lyotropic liquid crystals (Olmsted 2008) and numerous 
other systems. In the literature of non-Newtonian fluids (see, for a review, Olmsted 2008), the 
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shear-banding phenomenon has been explained as a constitutive instability from the linear sta- 
biUty analysis of appropriate constitutive models (Greco & Ball 1997; Wilson & Fielding 2005). 
The well-known "Hoffman-transition" in a colloidal suspension (banding/ordering of particles 
along the gradient direction) above the freezing density is accompanied by a sharp decrease in 
viscosity and has been explained in terms of a flow-instabiUty (Hoffman 1972). A very recent 
work (Caserta, Simeone & Guido 2008) on biphasic liquid-hquid systems showed that the shear- 
induced banding in such systems is tied to a lower viscosity, or, equivalently, a lower viscous 
dissipation. For both cases, the criterion of lower viscosity is similar to our criterion of a lower 
"dynamic" friction for the band-state in sheared granular fluid. It appears that the shear-induced 
banding in many complex fluids has a common theoretical description in terms of "constitutive" 
instability that leads to an ordered-state of a lower viscosity/friction. 
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